import os

import numpy as np
from scipy import signal
from scipy.signal import find_peaks, welch
from scipy.stats import iqr, variation
from vmdpy import VMD  # pip install vmdpy

# from vmdpy import VMD # 如果需要VMD，请安装此库并取消注释
def add(a, b):
    return a+b
def simple_entropy(p):
    """
    计算简单熵。
    """
    p = np.asarray(p)
    # 避免log(0)
    p = p[p > 0]
    E = -np.sum(p * np.log2(p))
    return E


def peak_check(angle, alpha, beta, fs, tempG1, wmod):
    """
    峰值检查函数。
    """
    Lcut = 4
    pos_p = 2  # Python中索引从0开始，原pos_p=3对应这里是2
    pos_z = 3  # 原pos_z=4对应这里是3

    # 使用welch替代pwelch
    fxx, Pxx_angle = welch(angle, window='hamming', fs=fs, nperseg=64, noverlap=32, nfft=64)
    Pxx_angle = Pxx_angle[Lcut:]

    # 查找峰值
    locs_welch, _ = find_peaks(Pxx_angle)
    if locs_welch.size > 0:
        id_angle = locs_welch[np.argmax(Pxx_angle[locs_welch])]
    else:
        id_angle = 0

    id_au = id_angle
    MPP = alpha * iqr(angle)
    # 确保fxx[id_au + Lcut - 1]不会越界
    if id_au + Lcut - 1 < len(fxx) and fxx[id_au + Lcut - 1] != 0:
        MPD = beta * 100 / fxx[id_au + Lcut - 1]
    else:
        MPD = 0.1  # 默认值，或者根据实际情况调整

    # 查找角度的峰值和谷值
    # MATLAB的findpeaks参数MinPeakProminence和MinPeakDistance
    # 在SciPy中分别对应prominence和distance
    pocks, _ = find_peaks(angle, prominence=MPP, distance=MPD, rel_height=0.5)
    locs  = angle[pocks]

    pocks_neg, _ = find_peaks(-angle, prominence=MPP, distance=MPD, rel_height=0.5)
    locs_neg = angle[pocks_neg]

    res = None  # 初始化 res 变量

    if tempG1 == 1:
        row1 = np.column_stack((pocks, angle[pocks], -np.ones(locs.shape), np.zeros(locs.shape)))
        row2 = np.column_stack((pocks_neg, angle[pocks_neg], np.ones(locs_neg.shape), np.zeros(locs_neg.shape)))

        res = np.vstack((row1, row2))

    elif tempG1 == 2:
        row1 = np.column_stack((pocks, angle[pocks], np.ones(locs.shape), np.zeros(locs.shape)))
        row2 = np.column_stack((pocks_neg, angle[pocks_neg], -np.ones(locs_neg.shape), np.zeros(locs_neg.shape)))

        res = np.vstack((row1, row2))

    res = res[res[:, 0].argsort()]  # 按第一列排序

    # 插入式均衡
    temp_flag = res[0, pos_p]
    temp_insertions = []
    temp_angle = angle

    for j in range(1, len(res)):
        if (temp_flag * res[j, pos_p] != -1):
            start_idx = int(res[j - 1, 0]) + 1
            end_idx = int(res[j, 0]) - 1

            if start_idx >= end_idx:  # 避免空区间
                continue

            cut_angle = temp_angle[start_idx:end_idx]

            if (temp_flag == -1 and tempG1 == 1) or (temp_flag == 1 and tempG1 == 2):  # 找谷值
                if cut_angle.size > 0:
                    Midx_rel = np.argmin(cut_angle)
                else:
                    continue
            elif (temp_flag == -1 and tempG1 == 2) or (temp_flag == 1 and tempG1 == 1):  # 找峰值
                if cut_angle.size > 0:
                    Midx_rel = np.argmax(cut_angle)
                else:
                    continue
            else:
                raise ValueError('tempFlag, tempG1错误')

            Midx = Midx_rel + start_idx
            temp_insertions.append([Midx, angle[Midx], -temp_flag, 0])
        else:
            temp_flag = res[j, pos_p]

    if temp_insertions:
        res = np.vstack((res, np.array(temp_insertions)))
        res = res[res[:, 0].argsort()]

    # 第二部分调整逻辑
    st = 1
    factor = 0.6
    # 确保wmod_temp和angle_temp是正确的形状
    wmod_temp = wmod / np.sqrt(np.mean(wmod ** 2))  # rms
    angle_temp = angle / np.std(angle, ddof=1)  # std
    i = 0
    for i in range(len(res) - 1):
        if tempG1 == 1:
            OCFlag = -1 * res[i, pos_p]
        else:
            OCFlag = 1 * res[i, pos_p]

        factor_t = factor * OCFlag

        end_idx = int(res[i + 1, 0])
        tempw_cut = -wmod_temp[int(st):end_idx]

        locs, _ = find_peaks(tempw_cut)
        pks = tempw_cut[locs]

        combined_metric = pks + factor_t * angle_temp[locs + st]

        if combined_metric.size > 0:
            id_max = np.argmax(combined_metric)
            Midx = locs[id_max] + st
            res[i, 0:2] = [Midx, angle[Midx]]
            st = int(res[i, 0])
        else:
            asdf = 1

    # 处理最后一个元素
    if res.size > 0:
        if tempG1 == 1:
            OCFlag = -1 * res[-1, pos_p]
        else:
            OCFlag = 1 * res[-1, pos_p]
        factor_t = factor * OCFlag
        temp_w_cut = -wmod_temp[st + 1:]
        if temp_w_cut.size > 0:
            locs, _ = find_peaks(temp_w_cut)
            pks = temp_w_cut[locs]
            combined_metric = pks + factor_t * angle_temp[locs + st]
            if combined_metric.size > 0:
                id_max = np.argmax(combined_metric)
                Midx = locs[id_max] + st + 1
                res[i+1, 0:2] = [Midx, angle[Midx]]

    # 删除重复或过近的索引
    delete_flag = 1
    while delete_flag == 1:
        delete_flag = 0
        delete_idx = []
        skip_count = 0
        for i in range(1, len(res)):
            if skip_count == 0:
                if res[i, 0] < res[i - 1, 0] + 6:
                    delete_idx.extend([i, i - 1])
                    skip_count = 2
                    delete_flag = 1
            else:
                skip_count -= 1
        res = np.delete(res, np.unique(delete_idx).astype(int), axis=0)

    # 检查交替性
    for i in range(1, len(res)):
        if res[i, pos_p] * res[i - 1, pos_p] == 1:  # 连续两个同类型峰值
            # 这部分在MATLAB代码中是assert(1==0)，表示不应该发生
            # 在Python中我们选择删除其中一个，或者抛出错误
            # 这里暂时选择抛出错误，如果需要更鲁棒的处理，可以修改
            raise ValueError("连续出现同类型峰值，请检查逻辑或数据")

    return res


def combine_with_wmod(angle, res, wmod, tempG1):
    pos_p = 2
    pos_z = 3
    threold = np.inf
    deleteIdx = []
    wmod = wmod / np.sqrt(np.mean(wmod ** 2))
    if angle.ndim > 1:
        angle = angle / np.std(angle, axis=0, ddof=1)
    else:
        angle = angle / np.std(angle, axis=0, ddof=1)

    st = 0  # 初始化 st
    for row_idx in range(len(res)):
        # 注意：res 中的索引（res[row_idx, 0]）在原始 Matlab 代码中可能是 1-based，
        # 所以在访问 wmod 时需要减 1
        if wmod[int(res[row_idx, 0]) - 1] < threold:
            st = row_idx  # 记录中断时的索引
            break
        else:
            deleteIdx.append(row_idx)

    # 如果 deleteIdx 不为空，从 res 中删除对应的行
    if deleteIdx:
        res = np.delete(res, deleteIdx, axis=0)
    assert st == 0

    idS = res[st, 0]
    skipFlag = 0
    w_in = []

    # 遍历 res
    for i in range(st + 1, len(res)):
        if skipFlag == 0:
            idE = res[i, 0]
            wmod_cut = wmod[int(idS):int(idE - 1)]
            inMax = 0  # 初始化 inMax
            if len(wmod_cut) > 0:
                inMax_relative = np.argmax(wmod_cut)
                inMax = inMax_relative + int(idS)
            else:
                inMax = int(idS)
            w_in.append([inMax, angle[int(inMax)], 0, res[i, pos_p]])
            idS = idE

    if w_in:
        w_in = np.array(w_in)
        assert len(w_in) == len(res) - 1

    if w_in is not None and w_in.shape[1] < res.shape[1]:
        padding = np.zeros((w_in.shape[0], res.shape[1] - w_in.shape[1]))
        w_in = np.concatenate((w_in, padding), axis=1)
    if res.shape[1] < 6:
        padding_cols_res = 6 - res.shape[1]
        padding_res = np.zeros((res.shape[0], padding_cols_res))
        res = np.concatenate((res, padding_res), axis=1)
    if w_in.shape[1] < 6:
        padding_cols_w = 6 - w_in.shape[1]
        padding_w = np.zeros((w_in.shape[0], padding_cols_w))
        w_in = np.concatenate((w_in, padding_w), axis=1)

    fin = np.concatenate((res, w_in), axis=0)
    fin = fin[fin[:, 0].argsort()]
    for i in range(len(fin)):
        if fin[i, pos_p] == -1 and fin[i, pos_z] == 0:
            fin[i, 5] = 1
        elif fin[i, pos_p] == 1 and fin[i, pos_z] == 0:
            fin[i, 5] = 3
        elif fin[i, pos_z] == 1 and fin[i, pos_p] == 0:
            fin[i, 5] = 2
        elif fin[i, pos_z] == -1 and fin[i, pos_p] == 0:
            fin[i, 5] = 4
        else:
            assert 1 == 0, "断言失败：不符合预期的条件"
    assert (np.sum((fin[:, 5] == 1) | (fin[:, 5] == 3))) == (np.sum((fin[:, 5] == 2) | (fin[:, 5] == 4)) + 1)
    for i in range(2, len(fin)):  # Python 索引从 0 开始，所以 i=3 对应 Python 的 i=2
        if (fin[i, pos_p] * fin[i - 2, pos_p] == -1 and fin[i, pos_z] * fin[i - 2, pos_z] == 0) or \
                (fin[i, pos_p] * fin[i - 2, pos_p] == 0 and fin[i, pos_z] * fin[i - 2, pos_z] == -1):
            pass
        else:
            assert 1 == 0, "断言失败：索引错位"
    deleteFlag = 1
    while deleteFlag == 1:
        deleteFlag = 0
        deleteIdx = []
        skipCount = 0

        # Matlab 的 j=2:length(fin(:,1))-2 对应 Python 的 range(1, len(fin) - 2)
        for j in range(1, len(fin) - 2):
            if skipCount == 0:
                # 注意 Matlab 的 fin(j-1,1) 对应 Python 的 fin[j-1, 0]
                # fin(j+2,1) 对应 fin[j+2, 0]
                idS = fin[j - 1, 0]
                idE = fin[j + 2, 0]

                if idE < idS + 2:
                    skipCount = 4
                    deleteFlag = 1
                    # Python 的 deleteIdx 存储的是行索引
                    deleteIdx.extend([j, j + 1, j + 2])

                    # Matlab 的 fin(j-1,4) 对应 Python 的 fin[j-1, 3] (pos_z)
                    if fin[j - 1, pos_z] == 0:  # 如果 j-1 为 wmod 峰值
                        # wmod(fin(j-1,1)) 对应 wmod[int(fin[j-1,0]) - 1]
                        # wmod(fin(j+3,1)) 对应 wmod[int(fin[j+3,0]) - 1]
                        if wmod[int(fin[j - 1, 0]) - 1] > wmod[int(fin[j + 3, 0]) - 1]:
                            deleteIdx.append(j + 3)
                        else:
                            deleteIdx.append(j - 1)
                    else:  # j-1 为 wmod 谷值
                        if j < 2:  # Matlab 的 j<3 对应 Python 的 j<2
                            deleteIdx.append(j - 1)
                        else:
                            # fin(j-1,1)-fin(j-3,1) 对应 fin[j-1,0]-fin[j-3,0]
                            # fin(j+3,1)-fin(j-1,1) 对应 fin[j+3,0]-fin[j-1,0]
                            if fin[j - 1, 0] - fin[j - 3, 0] > fin[j + 3, 0] - fin[j - 1, 0]:
                                deleteIdx.append(j + 3)
                            else:
                                deleteIdx.append(j - 1)
            else:
                skipCount = skipCount - 1

        # 删除重复的索引并排序，然后逆序删除以避免索引错位
        deleteIdx = sorted(list(set(deleteIdx)), reverse=True)
        for idx in deleteIdx:
            fin = np.delete(fin, idx, axis=0)

    # 最后的断言检查
    for i in range(2, len(fin)):
        if (fin[i, pos_p] * fin[i - 2, pos_p] == -1 and fin[i, pos_z] * fin[i - 2, pos_z] == 0) or \
                (fin[i, pos_p] * fin[i - 2, pos_p] == 0 and fin[i, pos_z] * fin[i - 2, pos_z] == -1):
            pass
        else:
            assert 1 == 0, "断言失败：最终索引错位"

    return fin


def arraylist_to_numpy(java_arraylist):
    size = java_arraylist.size()
    py_list = [java_arraylist.get(i) for i in range(size)]
    return np.array(py_list, dtype=np.float64)

def al_TH_PHI_w(w, fs, tempG1, X, Y, Z):
    print("Received w_T type:", type(w))
    print("Received X type:", type(X))

    # 将 Python 列表（从 Java 转换而来）转换为 NumPy 数组
    w = np.array(w)
    X = np.array(X)
    Y = np.array(Y)
    Z = np.array(Z)

    print("w_T shape:", w.shape)
    print("X shape:", X.shape)
    print("Y shape:", Y.shape)
    print("Z shape:", Z.shape)
    """
    主函数：计算特征。
    """
    alpha = 0.5
    beta = 0.4

    # 峰值存储索引 (Python中从0开始)
    pos_p = 2  # 原matlab中的3
    pos_z = 3  # 原matlab中的4

    wx = w[:, 0]
    wy = w[:, 1]
    wz = w[:, 2]
    wmod = np.sqrt(wx ** 2 + wy ** 2 + wz ** 2)

    res = peak_check(Z, alpha, beta, fs, tempG1, wmod)
    fin = combine_with_wmod(Z, res, wmod, tempG1)

    # 提取位置标签
    id_OP = fin[:, pos_p] == 1  # maxOpen
    id_CP = fin[:, pos_p] == -1  # maxClose
    id_maxO = fin[:, pos_z] == 1  # maxO
    id_maxC = fin[:, pos_z] == -1  # maxC

    idx_OP = fin[id_OP, 0].astype(int)
    idx_CP = fin[id_CP, 0].astype(int)
    idx_maxO = fin[id_maxO, 0].astype(int)
    idx_maxC = fin[id_maxC, 0].astype(int)

    # 切割出一段以闭合为开始，闭合为结束的拍打过程
    if idx_OP.size > 0 and idx_CP.size > 0:
        if idx_OP[0] < idx_CP[0]:
            idx_OP = idx_OP[1:]
            if idx_maxC.size > 0:  # 避免空数组操作
                idx_maxC = idx_maxC[1:]

        # assert((~any(idx_OP<idx_CP(1)))&&(~any(idx_maxC<idx_CP(1)))&&(~any(idx_maxO<idx_CP(1))));
        # Python中可以使用all()
        if idx_CP.size > 0:  # 确保idx_CP不为空
            assert (all(idx_OP >= idx_CP[0]) and all(idx_maxC >= idx_CP[0]) and all(idx_maxO >= idx_CP[0])), "初始索引检查失败"

        if idx_OP.size > 0 and idx_CP.size > 0:  # 再次检查空数组
            if idx_OP[-1] > idx_CP[-1]:
                idx_OP = idx_OP[:-1]
                if idx_maxO.size > 0:
                    idx_maxO = idx_maxO[:-1]
            assert (all(idx_OP <= idx_CP[-1]) and all(idx_maxC <= idx_CP[-1]) and all(
                idx_maxO <= idx_CP[-1])), "末尾索引检查失败"

    interval_array = np.diff(idx_CP)  # idx_CP(2:end)-idx_CP(1:end-1)
    if interval_array.size > 0:
        freq = 1. / interval_array
        freq_mean = np.mean(freq)
        SD_freq = np.std(freq, ddof=1)
    else:
        freq_mean = 0.
        SD_freq = 0.

    # 提取得到open过程和close过程
    open_process = []
    close_process = []

    if len(idx_CP) > 1 and len(idx_OP) > 0:  # 确保有足够的点进行配对
        # 确保idx_OP和idx_CP的长度匹配，或者进行适当的截断
        min_len = min(len(idx_CP) - 1, len(idx_OP))

        for i in range(min_len):
            open_process.append([idx_CP[i], idx_OP[i]])
            close_process.append([idx_OP[i], idx_CP[i + 1]])

    if open_process:  # 确保open_process不为空
        XYZ_X_range = np.abs(X[np.array([o[1] for o in open_process])] - X[np.array([o[0] for o in open_process])])
        XYZ_Y_range = np.abs(Y[np.array([o[1] for o in open_process])] - Y[np.array([o[0] for o in open_process])])
        XYZ_Z_range = np.abs(Z[np.array([o[1] for o in open_process])] - Z[np.array([o[0] for o in open_process])])
    else:
        XYZ_X_range = np.array([])
        XYZ_Y_range = np.array([])
        XYZ_Z_range = np.array([])

    assert len(open_process) == len(close_process), 'open,close长度失配'

    # 提取张开过程中的角速度
    OAV = []
    OAV_x, OAV_y, OAV_z = [], [], []
    OAV_first_cut, OAV_first_cut_x, OAV_first_cut_y, OAV_first_cut_z = [], [], [], []
    OAV_second_cut, OAV_second_cut_x, OAV_second_cut_y, OAV_second_cut_z = [], [], [], []

    N_open = len(open_process)
    N_open_cut = N_open // 2

    width_open = []
    for i, temp in enumerate(open_process):
        start, end = temp[0], temp[1]
        width_open.append(end - start)
        OAV.extend(wmod[start: end + 1])  # MATLAB的索引是包含end的
        OAV_x.extend(wx[start: end + 1])
        OAV_y.extend(wy[start: end + 1])
        OAV_z.extend(wz[start: end + 1])

        if i < N_open_cut:
            OAV_first_cut.extend(wmod[start: end + 1])
            OAV_first_cut_x.extend(wx[start: end + 1])
            OAV_first_cut_y.extend(wy[start: end + 1])
            OAV_first_cut_z.extend(wz[start: end + 1])
        elif i >= N_open - N_open_cut:  # MATLAB的N_open-N_open_cut+1 对应Python的 N_open-N_open_cut
            OAV_second_cut.extend(wmod[start: end + 1])
            OAV_second_cut_x.extend(wx[start: end + 1])
            OAV_second_cut_y.extend(wy[start: end + 1])
            OAV_second_cut_z.extend(wz[start: end + 1])

    OAV, OAV_x, OAV_y, OAV_z = np.array(OAV), np.array(OAV_x), np.array(OAV_y), np.array(OAV_z)
    OAV_first_cut, OAV_first_cut_x, OAV_first_cut_y, OAV_first_cut_z = np.array(OAV_first_cut), np.array(
        OAV_first_cut_x), np.array(OAV_first_cut_y), np.array(OAV_first_cut_z)
    OAV_second_cut, OAV_second_cut_x, OAV_second_cut_y, OAV_second_cut_z = np.array(OAV_second_cut), np.array(
        OAV_second_cut_x), np.array(OAV_second_cut_y), np.array(OAV_second_cut_z)

    # 提取闭合过程中的角速度
    CAV = []
    CAV_x, CAV_y, CAV_z = [], [], []
    CAV_first_cut, CAV_first_cut_x, CAV_first_cut_y, CAV_first_cut_z = [], [], [], []
    CAV_second_cut, CAV_second_cut_x, CAV_second_cut_y, CAV_second_cut_z = [], [], [], []

    N_close = len(close_process)
    N_close_cut = N_close // 2

    width_close = []
    for i, temp in enumerate(close_process):
        start, end = temp[0], temp[1]
        width_close.append(end - start)
        CAV.extend(wmod[start: end + 1])
        CAV_x.extend(wx[start: end + 1])
        CAV_y.extend(wy[start: end + 1])
        CAV_z.extend(wz[start: end + 1])

        if i < N_close_cut:
            CAV_first_cut.extend(wmod[start: end + 1])
            CAV_first_cut_x.extend(wx[start: end + 1])
            CAV_first_cut_y.extend(wy[start: end + 1])
            CAV_first_cut_z.extend(wz[start: end + 1])
        elif i >= N_close - N_close_cut:
            CAV_second_cut.extend(wmod[start: end + 1])
            CAV_second_cut_x.extend(wx[start: end + 1])
            CAV_second_cut_y.extend(wy[start: end + 1])
            CAV_second_cut_z.extend(wz[start: end + 1])

    CAV, CAV_x, CAV_y, CAV_z = np.array(CAV), np.array(CAV_x), np.array(CAV_y), np.array(CAV_z)
    CAV_first_cut, CAV_first_cut_x, CAV_first_cut_y, CAV_first_cut_z = np.array(CAV_first_cut), np.array(
        CAV_first_cut_x), np.array(CAV_first_cut_y), np.array(CAV_first_cut_z)
    CAV_second_cut, CAV_second_cut_x, CAV_second_cut_y, CAV_second_cut_z = np.array(CAV_second_cut), np.array(
        CAV_second_cut_x), np.array(CAV_second_cut_y), np.array(CAV_second_cut_z)

    maxCAV = wmod[idx_maxC]

    # 初始化
    fea = []
    nameTable = []

    # 迟缓相关特征计算
    # 速度的rms
    # 确保切片范围有效
    if idx_CP.size > 1:
        wmodCut = wmod[idx_CP[0]:idx_CP[-1] + 1]
        wxCut = wx[idx_CP[0]:idx_CP[-1] + 1]
        wyCut = wy[idx_CP[0]:idx_CP[-1] + 1]
        wzCut = wz[idx_CP[0]:idx_CP[-1] + 1]
    else:
        wmodCut, wxCut, wyCut, wzCut = np.array([]), np.array([]), np.array([]), np.array([])

    rms_wmod = np.sqrt(np.mean(wmodCut ** 2)) if wmodCut.size > 0 else 0.0
    rms_OAV = np.sqrt(np.mean(OAV ** 2)) if OAV.size > 0 else 0.0
    rms_CAV = np.sqrt(np.mean(CAV ** 2)) if CAV.size > 0 else 0.0

    # 拍打次数
    times = max(len(idx_OP), len(idx_CP)) - 1 if max(len(idx_OP), len(idx_CP)) > 0 else 0

    # 幅度计算
    MV_DR_XYZ_Y = np.mean(XYZ_Y_range) if XYZ_Y_range.size > 0 else 0.0

    # 幅度趋势分析
    N_range = len(XYZ_X_range)
    if N_range > 0:
        cut_length = N_range // 2
        first_cut = slice(0, cut_length)
        second_cut = slice(N_range - cut_length, N_range)

        CV_DR_XYZ_X = np.std(XYZ_X_range, ddof=1) / np.mean(XYZ_X_range)
        CV_DR_XYZ_Y = np.std(XYZ_Y_range, ddof=1) / np.mean(XYZ_Y_range)
        DE_DR_XYZ_X = np.sum(XYZ_X_range[second_cut] - XYZ_X_range[first_cut]) / cut_length if cut_length > 0 else 0.0
    else:
        CV_DR_XYZ_X = 0.0
        CV_DR_XYZ_Y = 0.0
        DE_DR_XYZ_X = 0.0

    CV_maxCAV = np.std(maxCAV, ddof=1) / np.mean(maxCAV)

    # 犹豫相关特征
    pks_wmod_neg, _ = find_peaks(-wmod)
    pks_wmod = wmod[pks_wmod_neg]  # 对应findpeaks(-wmod)

    # MATLAB的prctile函数，Python对应numpy.percentile
    if wmod.size > 0:
        # 计算 wmod 数组的第 25 个百分位数
        percentile_25 = np.percentile(wmod, 25)
        # 计算 pks 数组中小于 percentile_25 的元素数量
        count_less_than_percentile = np.sum(pks_wmod < percentile_25)
        # 计算最终结果 LMN_FAV
        LMN_FAV = count_less_than_percentile / times
    else:
        LMN_FAV = 0.0

    fea.extend([SD_freq, rms_wmod, rms_OAV, rms_CAV, times])
    nameTable.extend(['SD_Freq', 'RMS_FAV', 'RMS_OAV', 'RMS_CAV', 'times'])

    fea.extend([MV_DR_XYZ_Y, CV_DR_XYZ_X, CV_DR_XYZ_Y, DE_DR_XYZ_X, CV_maxCAV])
    nameTable.extend(['MV_DR_XYZ_Y', 'CV_DR_XYZ_X', 'CV_DR_XYZ_Y', 'DE_DR_XYZ_X', 'CV_maxCAV'])

    fea.extend([LMN_FAV])
    nameTable.extend(['LMN_FAV'])
    f, Cxy = signal.coherence(wx, wy, fs=1, window='hamming')

    peaks_xy_loc, _ = find_peaks(Cxy)
    peaks_xy = Cxy[peaks_xy_loc]
    Cxy_N = len(peaks_xy)

    # 能量熵特征
    if wxCut.size > 0:
        # 假设 wxCut, wyCut, wzCut 是列向量 numpy 数组
        # VMD 分解
        alpha = 1000  # 典型值，可调整
        tau = 0.  # noise-tolerance (no strict fidelity enforcement)
        K = 5  # 模态个数
        DC = 0  # no DC part imposed
        init = 1  # initialize omegas uniformly
        tol = 1e-7

        u, u_hat, omega = VMD(wxCut, alpha, tau, K, DC, init, tol)

        # 能量归一化
        swX = np.sum(u ** 2, axis=1)
        swX = swX / np.sum(swX)

        # 信息熵
        EwX = simple_entropy(swX)

        # PCA降维
        if wxCut.size >= 3:  # PCA需要至少3个样本
            # 这里的[wxCut,wyCut,wzCut]在MATLAB中是列向量组合成矩阵
            data_for_pca = np.column_stack((wxCut, wyCut, wzCut))
            # 手动PCA计算特征值（和MATLAB的lat类似）
            covariance_matrix = np.cov(data_for_pca, rowvar=False)  # rowvar=False表示列是变量
            eigenvalues = np.linalg.eigvalsh(covariance_matrix)  # 获取特征值
            eigenvalues = np.sort(eigenvalues)[::-1]  # 降序排列
            lat = eigenvalues / np.sum(eigenvalues)
            NPCV_FAV_1 = lat[0]
        else:
            NPCV_FAV_1 = 0.0  # 数据不足以进行PCA
    else:
        EwX = 0.0
        NPCV_FAV_1 = 0.0

    rms_wxCut = np.sqrt(np.mean(wxCut ** 2)) if wxCut.size > 0 else 0.0

    rate_fav_x = rms_wxCut / rms_wmod if rms_wmod != 0 else 0.0

    fea.extend([rate_fav_x, Cxy_N, EwX, NPCV_FAV_1])
    nameTable.extend(['Rate_Fav_x', 'Cxy_N', 'EwX_VMD', 'NPCV_FAV_1'])

    assert len(fea) == len(nameTable), "fea和nameTable长度未对齐"

    fea.extend([freq_mean])
    nameTable.extend(['MV_Freq'])

    return np.array(fea), nameTable



